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1.  INTRODUCTION 


Why  do  we  need  another  paper  on  the  estimation  of  point  and  mean 
free-air  gravity  anomalies  based  on  point  gravity  measurements?  Isn't 
this  subject  settled  once  and  for  all?  We  do  have  the  omnipotent  tool 
called  least-squares  collocation,  even  with  parameters. 

These  and  others  are  typical  questions  and  arguments  of  the  seven¬ 
ties,  when  the  geodetic  community  became  aware  of  and  excited  about  the 
existence  of  collocation.  Soma  enthusiastic  proponents  (not  its  de¬ 
signers,  mind  you)  advertised  it  as  the  unique  robot,  which  is  capable 
of  making  even  the  impossible  come  true;  but  unfortunately  it  doesn't. 
Meanwhile  the  enthusiasm  has  been  replaced  by  an  "unbiased"  recognition 
of  this  sound  and  powerful  tool;  practical  experience  has  shown  what 
could  have  been  anticipated:  a  smooth  and  reliable  output  requires  a 
smooth,  detailed  and  accurate  input  -  the  system's  response  is  data  - 
specific. 

One  of  the  many  applications  of  least-squares  collocation  is  the 
prediction  of  point  and  mean  gravity  anomalies  based  on  point  gravity 
anomalies  -  hardly  any  problem,  as  long  as  the  terrain  is  flat  within 
the  area  of  consideration.  Free-air  anomalies  can  be  processed  directly, 
no  data  reduction  seems  to  be  necessary.  If  we  approach  the  foothills 
or  even  mountainous  areas,  the  picture  changes  dramatically;  suddenly 
data  reduction  becomes  Indispensable,  collocation  results  based  on  un¬ 
reduced  quantities  become  practically  useless. 

This  report  aims  at  an  optimal  estimation  of  point  and  mean  anoma¬ 
lies  taking  into  account  the  concept  of  a  linear  correlation  between 
terrain-corrected  free-air  anomalies  and  topographic  height.  Different 
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estimation  procedures  are  compared;  particular  emphasis  is  put  on  the 
best  possible  estimation  of  the  Bouguer-factor  and  the  subsequent  estima¬ 
tion  of  Bouguer  as  well  as  free-air  point  and  mean  anomalies.  Least- 
squares  collocation  with  parameters  presents  itself  as  a  very  attractive 
and  powerful  tool  for  the  estimation  of  both  regression  parameters  and 
point  and/or  mean  anomalies.  An  explanation  for  the  regional  variation 
of  the  regression  parameters,  based  on  a  simplified  concept  of  isostatic 
compensation  is  presented  in  chapter  4. 

Particularly  in  mountainous  areas  the  free-air  anomalies  used  to 
be  reduced  for  the  effect  of  the  terrain,  which  can  easily  attain  values 
of  10-20  mgals  and  even  more.  !f  an  empirical  covariance  function  is 
estimated  from  unreduced  free-air  anomalies,  the  terrain  effect  causes 
the  variance  to  be  too  high  and  the  correlation  length  to  be  too  short. 
Since  the  variance  functions  as  a  scale  factor  for  the  prediction  error, 
we  see  that  the  quality  of  anomaly  prediction  is  worsened  if  unreduced 
free-air  anomalies  are  used.  In  addition,  the  correlation  length  controls 
the  quality  of  interpolation;  a  short  correlation  length  causes  a  large 
prediction  (interpolation)  error;  therefore,  the  prediction  accuracy  suf¬ 
fers  also  from  this  indirect  effect  (cf.  Siinkel ,  1981).  The  estimation 
of  the  terrain  effect  is  a  very  laborious  task.  It  is  shown  in  chapters 
5  and  6  to  depend  linear1/  on  the  terrain  variance  end  to  be  inverse 
proportional  with  respect  to  wavelength  and  correlation  length,  respec¬ 
tively.  The  total  variance  depends  on  the  power  spectrum;  if  the  high 
frequency  part  of  the  spectrum  has  much  power,  a  high  terrain  sampling 
rate  detailed  terrain  model)  is  required  in  order  to  estimate  the 
variance  with  sufficient  accuracy.  The  terrain  correlation  length  plays 
a  fundamental  role  in  the  estimation  of  the  mean  terrain  effect;  a  short 


correlation  length  requires  a  high  sampling  rate.  As  a  matter  of  fact, 
these  statistical  quantities  depend  strongly  on  the  terrain  in  consider 
ation;  consequently,  a  globally  valid  sampling  rate  cannot  be  given;  in 
dividual  circumstances  control  its  choice. 
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2.  MEAN  FREE-AIR  GRAVITY  ANOMALIES 

A  free-air  gravity  anomaly  is  defined  by 

Agp  *  9p  -  Yq  .  (2.1) 

with  gp  denoting  the  actual  gravity  at  the  point  P  located  on  the  sur¬ 
face  of  the  earth,  and  Yq  denoting  the  reference  gravity  at  a  correspond¬ 
ing  point  Q  located  on  the  telluroid  (Heiskanen  &  Moritz,  1967,  p.  293). 
By  Ag  we  denote  a  mean  value  of  the  free-air  gravity  anomaly, 

;  (2'2) 

Aop 

da  stands  for  the  element  of  solid  angle,  Ao  is  the  averaging  area  in 
consideration. 

Point  free-air  anomalies  are  known  to  oscillate  around  a  zero  aver¬ 
age  with  oscillation  frequencies  depending  strongly  on  the  topographical 
features,  they  comprise  regional  and  local  gravity  field  information; 
therefore,  a  point  anomaly  is  in  general  not  representative  for  a  large 
area;  because  of  its  local  character  and  its  strong  dependence  on  topo¬ 
graphy,  it  can  only  bo  poorly  predicted  if  the  topography  is  ignored 
(this  is  particularly  true  for  mountainous  areas). 

Averaging  free-air  anomalies  in  order  to  obtain  mean  anomalies  means 
essentially  averaging  the  local  features  of  free-air  anomalies;  mean 
anomalies  are  representative  for  the  averaging  area  in  consideration; 
they  offer  themselves  for  the  evaluation  of  various  integral  formulas. 

Appling  the  averaging  process  (2.2)  we  have  to  keep  in  mind  tr»at 
the  free-air  anomalies  refer  to  ground  level, 

Agp  a  Ag  (pp.Ap.h  {$p,Apj)  ;  (2.3) 


(h(<M)  denotes  the  topography.)  A  straight  mean  like  (2.2)  is  taken 
with  respect  to  the  horizontal  position  (<Mh  to  which  height  does  Agp 
refer?  Does  It  refer  to  a  mean  height? 

In  order  to  answer  this  question  we  introduce  a  mean  height  h  by 

*  AcT  ff  Vdcp-  •  (2-4) 

AOp 

Let  us  introduce  a  mean  anomaly  Ag*,  which  refers  to  the  mean  height 
h;  Ag*  is  a  mean  of  point  anomalies  Ag*,  which  refer  to  the  level  h  *  h 
and  are  obtained  by  an  anlytical  continuation  of  the  ground  level  anoma¬ 
lies, 

Ag*  -  Ag  -  ^  (h  -  h),  (2.5a) 

jf/Vv  •  '2-5b> 

Aop 

Denoting  the  average  (2.2)  by  M|  •  | ,  the  mean  anomaly  referred  to  the 
mean  height  becomes 

Ag*  -  Ag  +  h  M  j  ^  |  -M  j^hj  .  (2.6) 

It  is  instructive  to  consider  two  special  cases: 

a)  the  trivial  case  of  a  flat  topography  within  the  area  of  averaging, 
h  =  h;  it  is  obvious  that  the  last  two  terms  of  (2.6)  cancel  each 

other,  and 

Ag  *  Ag* 

follows.  The  straigt  mean  of  free-air  anomalies  refers  to  the  mean 
height. 

b)  the  vertical  gradient  3Ag/3h  is  constant  within  the  area  of  averag¬ 
ing,  3Ag/3h  ■  M  {3Ag/3h}  .  Again,  the  last  two  terms  of  (2.6)  can¬ 
cel  each  other,  and 
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follows.  The  straight  mean  refers  to  the  mean  height. 

The  last  case  Is  of  particular  Importance  because  It  assumes  a  linear 
relation  (•  exact  functional  dependence)  between  Ag  and  h.  In  reality, 
local  density  anomalies,  non-constant  Rouguer-anomalies  (in  the  averaging 
area),  nearby  topographic-  irregularities,  and  other  phenomena  account 
for  the  disturbance  of  an  exact  functional  dependence;  however,  in  gen¬ 
eral  we  observe  a  linear  correlation  (approximate  linear  relation)  be¬ 
tween  Ag  and  h,  which  becomes  even  more  ponounced  if  Ag  is  "cleaned"  from 
topographical  irregularities.  We  have  hardly  ever  a  sufficiently  dense 
gravity  material  which  would  allow  a  determination  of  the  vertical  gra¬ 
dient  cf  gravity;  therefore,  we  simply  have  to  assume  that 

H  (h-R)j 

is  zero,  in  other  words,  the  straight  mean  of  point  anomalies  Ag  is  inter¬ 
preted  as  Ag*  and  consequently,  the  reference  height  becomes  the  mean 
height.  In  areas  with  flat  topography  the  equality  holds  exactly. 


3.  PREDICTION  OF  ANOMALIES  BY  TREND  REMOVAL 


Adopting  the  concept  of  linear  correlation  between  free-air  anomaly 
and  topographic  height,  the  anomaly  can  be  represented  by 

Agp  -  a  +  bhp  +  sp  ,  (3.1) 

where  a  denotes  a  regional  constant  and  s  a  residual  anomaly;  s  comprises 
all  the  effects  which  make  Ag  locally  violate  the  linear  functional  re¬ 
lations;  the  average  of  s  is  assumed  to  be  zero, 

M  |s  f  *  0  .  (3.2) 

Then  the  average  of  the  mean  value  reduced  anomalies  Agr, 

Agp  =»  Agp  -  M  { Ag  | 

»  b  (hp  -  M  |h|  )  +  sp  ,  (3.3) 

vanishes,  M  {Agr}  3  0  . 

Introducing  a  reduced  height  hr, 

hp  =  hp  -  M  {  h  }  ,  (3.4) 

the  reduced  anomaly  Is  represented  as 

AgJ  *  bhj  +  sp  .  (3.3) 


If  s  is  to  be  independent  of  elevation,  it  follows  that  the  covariances 

W  MM 

cov  (Ag  ,  h  )  and  cov  (h  ,  h  )  are  proportional  to  each  other  for  all 
arguments  with  b  being  the  factor  of  proportionality  (Heiskanen  &  Moritz, 
1967,  p.  283  H.)  , 


b  =  cov  (Agr,  hr) 
cov  (hr  ,  hr) 


(3.5) 


Having  selected  a  best-fitting  b,  the  signal  s  can  be  predicted  at  any 
other  point  by  well-established  least-squares  prediction  methods.  A 
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free  air  anomaly  at  a  specific  point  Q  at  an  elevation  hg  is  obtained 
by 

Agg  *  M  { Ag  (  +  b  (hg  -  M  )  h  }  )  +  Sq  (3.6) 

or  shortly 

Agg  a  a  +  bhq  +  Sg  (3.6a)'' 

with  a  determined  through 

a  =  M  {  Ag  }  -  bM  {  h  }  .  (3.6b)' 

At  this  point  it  should  be  stressed  that  the  average  •  }  is  always 
derived  from  a  finite  sample  of  data  (free-air  anomalies,  topographic 
heights);  this  fact  will  not  cause  any  problems  as  long  as  the  terrain 
is  sufficiently  flat;  (the  average  over  a  perfectly  flat  terrain  is  ob¬ 
tained  by  a  single  data.)  However,  in  mountainous  areas  the  situation 
is  quite  different:  gravity  measurements  are  usually  performed  along 
main  leveling  lines  which,  in  turn,  almost  exclusively  coincide  with  main 
roads  running  through  valleys  rather  than  on  the  top  of  the  mountains. 
Therefore,  the  estimate  Mjh  }  will  tend  to  be  too  low,  and  since  the 
free-air  anomalies  are  linearly  correlated  with  height,  the  estimate 
M  {  Ag  }  tends  to  be  foo  low  as  well.  The  direct  effect  on  the  point 
anomaly  prediction  would  be  negligible  provided  the  gravity  gradient  b  is 
sufficiently  well  determined.  However,  the  estimation  of  b  is  poor  if 
the  height  range  of  the  data  is  narrow.  Therefore,  sampling  gravity  and 
corresponding  height  in  valleys  only  will  make  point  predictions  very 
inaccurate.  A  homogeneous  sampling  is  therefore  strongly  recommended. 

A  rrv-an  anomaly  Ag  can  be  derived  immediately  from  (3.6), 

Ag  =  M  {  Ag  }  +  b  (h  -  M{  h  }  ), 


(3.7) 
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provided  the  average  of  the  signal  s  vanishes.  Again  we  observe  that 
the  uncertainty  of  b  enters  directly  in  the  uncertainty  of  the  mean  value. 
Practical  results  (Uotila,  1967  a,b;  SUnkel  und  Malits,  1981)  indicate 
that  b  can  be  estimated  with  an  error  ranging  from  ±10' 3  to  ±10~2  mgal/m; 
therefore,  a  difference  Ah  =  h  -  M{  hf  of  a  couple  of  hundred  meters 
(a  case  quite  likely  In  mountainous  areas)  can  contribute  to  the  error 
budget  of  Ag  by  a  couple  of  mgals. 

As  far  as  the  estimation  of  a  and  b  is  concerned,  either  a  classical 
least-squares  fit  (Uotila,  196/  a,b)  or  a  more  demanding  least-squares 
collocation  solution  can  be  envisioned.  The  least-squares  estimation 
of  the  parameters  a  and  b  (trend  model  parameters)  is  performed  by  con¬ 
sidering  s  as  random  noise.  Denoting  the  parameter  vector  by  X  and  the 
design  matrix  by  A, 


_  r  i,  i . i 

(  =  (a,  b)  ,  At  *  h  h  h 

n1,n2,. ..,  nn 


the  best  estimate  of  X  is  obtained  by 


X  =  (ATA)_1ArAg; 


(3.8) 


its  elements  can  easily  be  shown  to  equal 

n  r 

1  hi  Ag, 

e  =  iii -  ,  a  5  M  J  ig  }  -bM(hr|  .  (3.8)' 

"  y*  2 

2(h,) 

i  =  l 

The  least-squares  collocation  solution  differs  from  the  simple  least- 
squares  solution  insofar  as  it  also  takes  into  account  the  statistical 
behavior  of  the  signal  s  (which  is  for  sure  not  simply  random  noise), 
expressed  in  terms  of  its  covariance  function.  The  best  estimate  of  the 
parameter  vector  X  is  now  given  by 
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X  *  (AtC_1A}  1  ATC"1Ag  (3.9) 

where  C  denotes  the  covariance  matrix  of  the  vector  of  signals.  There¬ 
fore,  the  two  solutions  will  differ,  if  the  individw1  r  are  cor¬ 
related- (C  has  non-vanishing  off-diagonal  elements.)  It  need  not  to 
be  mentioned  that  the  collocation  solution  (3.9)  is  by  far  more  expen¬ 
sive  than  the  least-squares  solution  (3.8),  because  it  requires  a)  an 
estimation  of  an  empirical  covariance  function  of  the  signal  r ,  b)  the 
fit  of  an  appropriate  model  to  the  empirical  covariance  function,  c)  the 
calculation  of  [n  (n-1)]  /2  covariances,  and  d)  the  inversion  of  the  co- 
variance  matrix.  This  is  the  price  we  have  to  pay  for  an  optimal  esti¬ 
mate  obtainable  from  the  available  data  set.  The  price  will  be  defi¬ 
nitely  too  high  if  the  correlation  between  the  signals  is  very  weak 
(almost  diagonal  covariance  matrix),  it  will  be  a  good  investment  if 
the  correlations  are  strong:  practical  studies  have  shown  (S'unkel  and 
Malits,  1981)  that  the  variance  of  the  signal  can  be  quite  considerable 
(easily  reach  100  mgal  and  more);  this  variance  enters  unreduced  into 
the  error  estimate  of  the  point  anomaly  if  the  simple  least-squares  con¬ 
cept  is  applied,  in  the  least-squares  collocation  method  the  predicted 
signal  is  obtained  from  the  centered  da*::  by  the  well-known  relation 

§Q  =  Cj  C"1  (Ag  -  AX);  (3.10) 

the  error  covariance  matrix  Ess  of  the  predicted  signal  vector  s  is  given 
b.y 

Ess  =  Css  -  C^C-1  [i  -  A(ATC‘1A)“1ATC'1]  Cs  ,  (3.11) 

the  error  covariance  matrix  of  the  predicted  parameter  vector  by 
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Exx  =  (A'C'^A)”1  (3.12) 

(Moritz,  1900,  p.  128).  Equation  (3.11)  can  be  split  up  into  three 
terms, 

^ss  =  css  "  CSC  Cs  +  (A  C  CS)EXX  (A  C  Ci)  ; 

the  first  term  represents  the  a  priori  error  (covariance  matrix)  of  the 

estimated  signal  vector,  the  second  term  represents  the  accuracy  gain 
due  to  data,  and  the  tiird  term  the  contribution  of  the  parameter  inac¬ 
curacies  to  the  error  of  the  estimated  signal.  As  a  general  rule  it 
can  be  said  that  a  strong  signal  variance,  a  small  correlation  length 

compared  to  the  mean  mutual  distance  between  data,  and  a  vague  linear 

correlation  between  Ag  and  h  will  be  responsible  for  poor  signal  ac¬ 
curacies;  a  low  variance,  long  correlation  length,  and  a  strong  linear 
correlation  will  keep  the  prediction  error  low. 

Any  other  free-air  point  gravity  anomaly  can  be  estimated  through 

AQ0  *  V  +  §o  (3.13) 

where  A  denotes  the  design  vector  corresponding  to  the  prediction 
point  Q, 

a;  -  (1,  h0)  .  (3.14) 

The  mean  square  estimation  error  is  obtained  by 

m^gQ  55  AQEXXAQ  +  AqExs  +  ExsAq  +  Ess  (3.15a/ 

with  the  error  cross  covariance  matrix 


Exs  -  -(A,C"1A)"lA,C"1Cs 


(3.15b) 


msmgmssgsBim 
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Equation  (3,13)  enables  us  to  estimate  the  mean  free-air  anomaly  through 
Ag  =  A0X  +  s0  ,  (3.16) 

where  A0  denotes  the  mean  of  all  vectors  (3.14)  which  is  obviously 
A0  «  (1,  fi)  . 

The  estimated  mean  anomaly  refers  to  the  mean  height;  if  the  mean  of 
the  signal  vanishes,  we  obtain 

»  3  +  6  h  (3.17) 

as  best  linear  unbiased  estimate  of  the  mean  gravity  anomaly.  (Note 
that  it  makes  no  difference  if  we  take  the  mean  of  the  surface  anoma¬ 
lies  or  the  mean  of  the  anomalies  at  mean  height;  this  is  because  the 
signal  was  supposed  to  be  uncorrelated  with  height.)  Formally  (3.17) 
is  valid  for  both,  the  adjustment  solution  and  the  least-squares  col¬ 
location  solution;  however,  the  way  the  estimates  i  and  b  are  obtained 
is  different. 

The  collocation  solution  is  based  on  the  covariance  function  of  the 
signal  s;  therefore,  it  would  oe  quite  interesting  to  know  which  kind 
of  information  is  contained  in  s.  According  to  its  definition  (3.1), 
s  is  a  kind  of  mean  value  reduced  Bouguer  anomaly  with  the  parameter  a 
as  the  (constant)  mean  Bouguer  anomaly  in  the  area  of  consideration  and 
b  the  Bouguer  factor  which  can  be  related  to  the  mean  density  p  by  b  = 

271  Gp  (G  is  the  gravitational  constant).  Since  no  terrain  correction  has 
been  taken  into  account  in  our  model,  the  signal  s  will  comprise  es¬ 
sentially  3  kinds  of  information: 


mam* 
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a)  the  variation  of  the  Bouguer  anomaly  within  the  area  of  con¬ 
sideration, 

b)  the  terrain  effect, 

c)  the  effect  of  local  and  regional  density  anomalies. 

Terrain  corrected  Bouguer  anomalies  Ag#  are  known  to  be  very  smooth  and 
correlated  to  a  "mean"  height  h  such  that  in  average 

Ag#  3  -  100  •  ft  [  km  ]  mgal  (3.18) 

(Heiskanen  &  Moritz,  1967,  p.  328);  the  definition  of  "mean"  is  strongly 
linked  to  the  concept  of  Isostasy  which  will  be  discussed  in  the  next 
chapter.  If  the  terrain  is  rough,  the  behavior  of  the  terrain  correc¬ 
tion  will  be  similarly  rough  and  high-frequent.  Since  no  terrain  cor¬ 
rection  has  been  applied,  the  signal  s  will  have  a  long-wavelength  Bouguer 
anomaly  characteristic  superimposed  by  a  short-wavelength  terrain  effect 
characteristic;  density  anomalies  will  probably  cover  the  whole  spectrum 
as  far  as  its  effect  on  the  signal  is  concerned  however,  its  power  is 
considered  to  be  rather  small  compared  to  the  Bouguer  and  terrain  effect. 

In  moderately  rough  areas  the  terrain  effect  will  be  small  and  the 
signal  is  essentially  a  Bouguer  anomaly.  It  is  true  that  Bouguer  anoma¬ 
lies  are  smooth  (ci$.  (3.18)),  but  still  they  are  hardly  ever  constant 
in  an  area  of  say  1°  x  1°;  its  variation  enters  fully  into  the  signal  s. 
What  is  the  impact  of  the  relation  (3.18)  on  a  determination  of  the  para¬ 
meter  b  by  least-squares  adjustment?  Very  simple,  the  value  of  b  will 
tend  to  decrease  with  increasing  area  of  consideration.  Why?  As  we 
saw  earlier,  s  is  essentially  a  Bouguer  anomaly  which  behaves  accord- 
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^  P  . 

$*  Ag,  +  const,  a-O.l'h  [m  j  +  const., 

and  consequently,  the  free-air  anomaly  is  approximately  given  by 

Agf »  bhf  -  O.lh,  +  corst., 

and  with  homogenous  density  p  =  2.67  g  cm”3, 

Ag,=  0.1(h  -  h),  +  const.  .  (3.19) 

Consider  a  region  of  1°  x  1°,  subdivide  it  into  9  subregions  of  size 
20*  x  20' ,  and  assume  that  the  mean  height  h  is  constant  within  each 
subregion,  but  changes  from  region  to  region  (a  sufficiently  justified 
assumption  as  shown  in  the  next  chapter);  assume  furthermore  that  dense 
gravity  material  is  available  in  the  whole  region.  Linder  these  assump¬ 
tions  we  plot  the  free-air  anomalies  versus  height  for  each  sub-region 
and  should  expect  a  correlation  behavior  which  shows  a  parallel  shift 
from  one  sub-region  to  the  next  as  illustrated  in  Figure  3.1  below. 


Fig.  3.1  Correlation  model  between  free-air 

anomaly  and  topographic  height.  j 

i 

i 

s 

I 


I 
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(note  that  the  parallel  shift  is  due  to  the  different  mean  heights.)  If 
we  perform  a  least-squares  adjustment  solution  for  the  estimation  of  the 
parameter  b  for  the  whole  region,  we  will  get  a  too  small  value  b  which 
is  Indicated  by  the  direction  of  the  boldface  arrow  In  Fig.  3.1.  The 
reasoning  is  very  simple:  the  parallel  shift  is  essentially  the  signal 
s;  in  the  adjustment  procedure  this  signal  is  treated  as  random  noise; 
with  other  words,  the  adjustment  solution  is  blind  with  respect  to  hori¬ 
zontal  position,  it  simply  makes  a  mishmash  of  Ag  with  respect  to  h  and 

t 

not  of  its  gradient  as  it  should. 

In  contrast  to  the  adjustment  solution,  the  least-squares  colloca¬ 
tion  solution  takes  care  of  the  signal  very  carefully;  it  uses  all  the 
information  contained  in  s,  provided  in  terms  of  covariances,  in  order 
to  estimate  a)  a  common  gradient  b  and  b)  the  signal  field  as  such  (equ. 
(3.9),  (3.10)).  Therefore,  a  regional  collocation  solution  will  signi¬ 
ficantly  improve  the  estimation  of  b  as  compared  to  a  regional  adjustment 
solution  whenever  a)  the  region  is  large,  and  b)  the  terrain  is  not  flat. 

It  should  be  mentioned  that  Uotila  ( 1967a, b)  performed  extensive 
numerical  studies  in  order  to  find  an  optimal  procedure  for  the  estima¬ 
tion  of  a  regionally  valid  parameter  b.  He  finally  came  to  the  conclu¬ 
sion  that  a  resonable  estimate  can,  in  general,  only  be  achieved  if  the 
region  Is  subdivided  into  smaller  subblocks;  for  each  subblock  a  local 
parameter  b  should  be  estimated  by  least-squares  adjustment,  the  regional 
value  is  obtained  as  an  appropriate  average  of  all  local  b  -  values. 

His  conclusions  are  in  favorable  agreement  with  ours;  here  a  simple 
explanation  of  this  phenomenon  has  been  provided.  A  m; thematically  sound 
reasoning  will  be  attempted  in  chapter  4. 
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Let  us  turn  back  to  the  signal  and  its  statistical  behavior.  The 
collocation  estimation  of  the  parameter  vector  (3.9)  requires  the  know¬ 
ledge  of  the  signal's  covariance  function  which  we  do  net  necessarily  have 
at  hand;  however,  if  a  sufficiently  large  numoer  of  data  is  available 
within  the  considered  region,  we  can  proceed  iteratively:  we  chose  the 
standard  value  b  *  0.112  as  starting  parameter,  obtain  the  signal  at 
the  data  points  by 

sf  ■  Ag,  -  M  |Ag|  -  b0  (h,  -  M{h}  ),  (3.20) 

calculate  an  empirical  covariance  function,  and  fit  an  appropriate  model 
covariance  function  which  is  used  in  the  estimation  of  the  parameter 
vector.  If  necessary  (hardly  ever  it  is)  we  can  determine  better  esti¬ 
mates  of  the  signal  and  so  forth.  By  far  the  most  expensive  part  in  the 
least-squares  collocation  solution  remains  the  calculation  of  the  indivi¬ 
dual  covariances  and  the  inversion  of  the  covariance  matrix. 

The  least-squares  collocation  solution  presented  here  was  based  on 
the  very  essential  and  restricting  assumption  that  free-air  anomalies 
and  elevations  are  linearly  correlated;  in  other  words,  we  have  assumed 
that  the  covariance  functions  cov(Agr,  hr)  and  cov(hr,  hr)  are  propor¬ 
tional  for  all  arguments,  a  condition  which  was  imposed  in  order  to  render 
the  signal  s  independent  of  elevation.  However,  if  this  condition  is 
intolerably  violated,  a  more  general  approach  has  to  be  aimed  at.  In  a 
very  early  paper  Moritz  (1963)  has  laid  down  the  basic  formulas  which 
express  the  optimally  predicted  centered  free-air  anomaly  in  terms  of 
a  linear  combination  of  centered  free-air  anomalies  and  centered  eleva¬ 
tions.  The  predicted  quantity  is  given  by  the  familiar  expression 


A9:  •  eje-1  {  . 


(3.21) 
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In  this  context  i  denotes  the  vector 

it  s  (Ag^  ,  Ag2  i  . .  • ,  Agn  »  h,  *  h2  ,  •  • • »  hR  »  hf )  > 

note  that  the  elevation  hj  of  the  prediction  point  is  an  element  of  the 
"data"  vector.  This  peculiar  case  deserves  some  attention:  the  free- 
air  anomaly  is  a  ^unction  of  horizontal  and  vertical  position;  its  re¬ 
striction  to  the  surface  of  the  earth  Is  characterized  by  (2.3).  A 
linear  predictor  has  to  represent  Ag  in  terms  of 

Ag^  *  o£  Agr  +  B*hr  +  yph*  (3.22) 

with  coefficients  aj,  and  y,  independent  of  elevation.  Since  hr  is 
correlated  with  Agr,  the  minimization  of  the  prediction  error  leads  to 
a  linear  system  which  consists  of  covariances  between  data  (usual  case) 
and  covariances  between  the  elevation  at  the  prediction  point  and  all 
data.  Therefore,  we  obtain  a  covariance  matrix  which  can  be  partitioned 
into  4  blocks. 


with  the  prediction-point-independent  (data,  data)  -  covariance  matrix  C2 
anu  the  variance  of  centered  heights  C3,  and  the  prediction-point-dependent 
block  C2  which  is  a  vector  of  covariances  between  the  height  at  the  pre¬ 
diction  point  and  all  data.  When  we  are  talking  about  "data",  we  have 
the  set 

(Agi  ,  Ag2  ,  . ..,  Ag„  ;  hj  ,  h2  ,  ...»  h„) 


in  mird. 
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Here  we  see  already  the  drawback  of  this  general  and  optimal  pre¬ 
diction  method: 

A)  We  need  to  know 

a)  the  autocovariance  function  of  free-air  anomalies  cov(Agr,Agr) , 

b)  the  autocovariance  function  of  the  topography  cov  (hr,hr), 

c)  the  crosscovariance  function  of  free-air  anomaly  and  topography; 

B)  Each  prediction  requires  the  calculation  and  inversion  of  a  co- 
variance  matrix  which,  in  contrast  to  usual  least-squares  predic¬ 
tion  problems,  changes  with  the  horizontal  position  of  the  predic¬ 
tion  point;  (the  covariance  matrix  is  invariant  with  respect  to  the 
vertical  position  of  the  prediction  point  only,  but  not  with  respect 
to  its  horizontal  position.) 

As  a  matter  of  fact,  this  peculiar  property  of  the  covariance  matrix 
makes  predictions  more  expensive,  however,  not  as  much  as  one  would  expect 
from  a  first  glance,  for  the  following  reason:  the  block  Cx  is  invariant 
with  respect  to  the  location  of  the  prediction  point  and,  therefore,  Ci 
has  to  be  inverted  only  once  and  for  all.  The  remaining  parts  of  the 
inverse  covariance  matrix  can  be  obtained  by  a  simple  block-partitioning 
( C(5 .  Faddeeva,  1959,  §14): 


r  c i  c2Vl  r  Bi  b2‘ 

=*■  ■  I  ■ 

.  C2  C3  j  1.  B;  B3. 

with  B3  -  1/(C3  -  C2TClxC2)"1 
B2  *-  Ci  C2Bj 


(3.23) 


Bi  =  Cll+  B282/Bj 
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Assuming  n  data  given  (1  data  consists  of  Agr  and  the  corresponding  hr), 

Ci  has  dimension  2n  x  2n,  C2  has  dimension  2n  x  1,  and  Cj  is  a  constant. 
Observing  (3.23)  we  conclude  that  on  the  order  of  8n2  basic  operations 
(multiplication  +  addition)  are  required  for  the  calculation  of  the  full 
inverse  covariance  matrix  provided  that  Cl1  is  already  available.  In 
viaw  of  the  fact  that  the  calculation  of  Cj1  requires  on  the  order  of  8n3 
basic  operations,  we  conclude  that  the  prediction  of  n  anomalies  is  just 
twice  as  expensive  as  the  calculation  of  Cl1  (which  Is  approximately  equal 
to  the  prediction  of  a  single  point  anomaly).  Therefore,  the  dependence 
of  the  covariance  matrix  on  the  horizontal  location  of  the  prediction 
point  will  not  blow  up  the  computation  time  too  much  and  cannot  be  con¬ 
sidered  a  severe  limitation.  There  are  much  stronger  arguments  which  do 
not  speak  in  favor  of  this  "optimal"  solution. 

Let  us  compare  the  collocation  solution,  which  assumes  linear  cor¬ 
relation  between  Ag  and  h,  with  the  collocation  solution  based  on  a 
rather  arbitrary  correlation  behavior.  The  solution  (3.10)  requires 
the  estimation  of  a  single  covariance  function,  its  fit  by  an  appropriate 
model,  the  calculation  of  about  n2/2  covariances,  and  the  inversion  of  a 
n  x  n  symmetric  and  positive  definite  matrix  if  n  gravity  anomalies  are 
processed.  The  general  method  requires  the  estimation  of  3  times  more 
covariance  functions,  their  fit  by  appropriate  models,  the  calculation 
of  about  4  times  more  covariances  and  8  times  more  basic  operations  for 
the  calculation  of  the  Inverse  covariance  matrix.  Therefore,  it  seems 
to  be  on  the  order  of  5  times  more  expensive  than  the  collocation 
solution  based  on  the  linear  correlation  model.  Particularly  trouble¬ 
some  seems  to  be  the  estimation  of  cov(Agr,hr)  and  cov(hr,hr)  because 
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of  the  generally  lacking  data  density*  particularly  in  mountainous  areas, 
where  there  would  be  a  real  need  for.  For  these  reasons,  it  Is  rather 
questionable  If  the  optimal  solution  for  the  general  case  differs  signi¬ 
ficantly  from  the  linear  correlation  solution  (if  a  linear  correlation 
exists,  both  solutions  are  Identical),  and  if  a  small  improvement  of  the 
solution  justifies  the  very  high  price  to  be  paid. 

The  following  Table  summarizes  rms  -  prediction  errors  of  mean  free- 
air  anomalies,  depending  on  the  data  density  and  the  correlation  length 
£  of  the  covariance  function;  a  variance  of  100  mgal2  has  been  used.  The 
data  are  assumed  to  he  regularly  distributed  error-free  point  gravity 
anomalies.  A  data  density  of  N  means  N  data/ blocks. 

The  figures  In  Table  3.1  do  not  include  the  error  introduced  by  the 
inaccuracies  of  the  trend  model  parameters  a  and  b;  they  can  contribute 
to  the  total  error  budget  up  to  2  -  3  mgal  (Silnkel  and  Mai  its,  1981); 
a  typical  error  estimate  of  a  is  ±  1.5  mgal,  a  typical  error  estimate 
of  b  is  ±0.002  mgal/m.  Note  that  the  variance  of  100  mgal2  has  been 
chosen  rather  arbitarlly;  the  figures  can  easily  be  scaled  by  the  proper 
variance. 
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Table  3.1  Mean  free-air  anomaly  rms  -  prediction  errors  (mgal),  depend¬ 
ing  on  the  data  density  N  (*  number  of  data/block)  and  on  the 
correlation  length  £;  common  variance:  100  mgal2. 


The  results  summarized  in  Table  3.1  are  graphically  represented  in  Fig. 
3.2  a-e.  The  contours  are  lines  of  constant  mean  free-air  anomaly  predic 
tion  error  dependent  upon  the  correlation  length  £  (horizontal  axis)  and 
the  data  density  (vertical  axis).  Note  that  these  estimates  refer  to  the 
ideal  situation  of  regularly  distributed  and  error-free  data  having  a 
variance  of  100  mgal2; 
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I, 

r  Fig.  3.2  a-e  Lines  of  constant  rms  mean  anomaly  prediction 

I  error,  depending  on  the  correlation  length  and 

i  the  data  density;  variance:  100  mgal2 


i  As  a  matter  of  fact  the  estimates  obtained  here  have  to  be  scaled  accord- 

1  ing  to  the  individual  variance.  The  variance  can  considerably  differ 

K' 

r 

^  between  various  areas;  this  is  why,  as  a  kind  of  normalizing  factor,  a 

\  common  100  mgal2  variance  has  been  assumed.  E.g.  the  mean  residual  var- 

I  iance  for  a  1°  x  1°  area  is  of  the  order  of  900  mgal2  (total  variance 

i 

l  minus  variance  of  1°  x  1°  mean  anomalies);  therefore,  the  values  of  Fig. 

r 

I  3.2 e  should  (in  average)  be  multiplied  by  a  factor  3.  All  values  refer 

to  the  covariance  model  of  Hirvonen  C(s)  =  C0/[l+(s/£)2] ;  for  the 

l  a  2  -  2 

i  Gaussian  covariance  model  C(s)  *  Coe  we  obtained  estimates  which 

I  are  about  50%  lower  for  S'  x  S'  mean  values  and  15%  lower  for  1°  x  1°; 

these  lower  estimates  are  due  to  the  stronger  correlation  of  the 
\  Gaussian  model,  compared  to  the  Hirvonen  model,  for  small  distances;  the 

\  covariance  function's  behavior  for  small  distances  is  controlled  by  its 

I 

|  curvature  parameter  at  the  origin;  therefore,  highly  reliable  prediction 


and  prediction  error  estimates  should  be  based  on  a  covariance  model  which 
resembles  all  three  essential  parameters,  the  variance,  the  correlation 
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4.  The  b  =  const.  -  PROBLEM  and  ISOSTASY 

In  chapter  3  basically  three  methods  for  the  prediction  of  free- 
air  anomalies  have  been  presented.  Two  of  them  assumed  a  linear  cor¬ 
relation  between  free-air  anomalies  and  elevations,  represented  by  the 
correlation  coefficient  b.  In  least-squares  adjustment  determinations 
of  b  (equ.  (3.8)')  ,  it  has  been  observed  (Uotila,  1967  a,b)  that  |b| 
tends  to  decrease  if  the  area,  for  which  it  is  considered  constant, 
increases.  A  simple  explanation  of  this  phenomenon  has  been  provided; 
it  was  based  on  the  assumption  that  the  value  of  the  Bouguer  anomaly  is 
approximately  proportional  to  a  "mean"  height.  It  has  been  anticipated 
that  the  way  of  taking  the  mean  of  the  topography  is  closely  linked  to 
the  concept  of  Isostasy.  In  this  chapter  we  make  the  attempt  to  obtain 
a  mathematical  relation  between  the  area  size  (of  b  considered  constant), 
statistical  characteristics  of  the  topography,  and  the  error  to  be  ex¬ 
pected  in  b  due  to  the  b  *  const,  assumption. 

Moritz  (1969)  has  shown  by  means  cf  a  simplified  model  of  isostasy, 
that  the  linear  correlation  between  free-air  anomalies  and  topographic 
elevations  can  be  explained  in  a  very  simple  way.  Representing  the 
compensated  masses  by  a  surface  layer  at  the  depth  D  below  sea  level, 
he  finally  arrives,  after  some  minor  neglections,  at  a  relation  which 
expresses  the  free-air  anomaly  In  terms  of  the  corresponding  point 
height,  the  corresponding  mean  height,  and  the  topographic  correction. 


Agp  =  2ttGo  (h,  -  hf)  -  C, 


(4.1) 
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h,  ...  elevation  of  the  point 

ft,  ...  mean  elevation  corresponding  to  P, 

C,  ...  topographical  correction, 

p  ...  density, 

G  ...  gravitational  constant. 

(Note  that  for  this  model  the  Bougucr  anomalies  are  given  by  ~?TrGpftp , 
the  isostatic  anomalies  vanish.) 

The  mean  height  ft,  is  represented  in  terms  of  the  output  of  a  linear 
system  witn  input  h, 

ft(p)  a  YT  j j  £4^)  da  (Q)  (4.2) 


with  lt  denoting  the  distance  between  the  point  P0  (located  at  sea  level, 
orthogonal  projection  of  the  surface  point  P)  and  Q  (located  on  the 
compensation  surface  at  depth  D  below  sea  level);  R  denotes  the  mean 
radius  of  the  earth.  In  the  spectral  domain  the  mean  height  spectrum 
ft**,  is  related  to  the  point  height  spectrum  hww  through 

h,m  ®  ,  (4.3) 


where  k,  is  the  n'th  degree  eigenvalue  of  the  integral  kernel  K  of  (4.2), 


K{P.Q):  a  • 

According  to  the  Funck-Hecke  formula  (Miiller,  1966,  p.  20), 
K.  =  2tt  J K(t)P.  (t)dt; 


(4.4) 


(4.5) 
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(?„  is  the  Legendre  polynomial  of  degree  n.)  Introducing  (4.4)  in  (4.5), 
k.  is  expressed  by 


-l 


as  the  distance  between  P0  and  Q  is  given  by 
it  -  [r2  +  (R-D)2  -  2R  (R-D)  cos^]**, 
and  1 fl\  by 

11 

jt  *  p-  (1  +  <*2  -  2  a  t)  (4.6a) 

with  a  :  !l-|,t:  *  cos^  .  (4.6b) 

The  expression  (4.6a)  can  be  represented  in  terms  of  a  series  of  Legendre 
polynomials  (Heiskanen  &  Moritz,  1967,  p.  35), 

j r  *  "RT(  1-a2 )  X  ^2n  +  ^  P"  (cos^  •  (4--') 


Now  it  is  fairly  easy  to  derive  the  eigenvalues  ;  observing  the  ortho 
gonality  relation  of  Legendre  polynomials  expressed  through 

i 

j  ”,  (t)  pm  (t)  dt  =  smm 

(6b.  denotes  the  Kronecker  symbol),  we  obtain  with  (4.5)  and  (4.7) 


< 


2ot" 

n^r  ; 


(4.8) 
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expressing  a  by  (4.6b)  we  obtain 


In  order  to  better  understand  the  impact  of  these  eigenvalues,  let  us 
consider  two  extreme  cases: 

a)  D  «  0  (compensation  level  coincides  with  sea  level):  in  this  case 
both  the  isostatic  and  the  Bouguer  anomalies  should  be  identically 
zero  according  to  our  model  (4.1).  In  other  words,  fr  i  h  -  no 
smoothing  is  envoi ved.  In  terms  of  the  spectrum  this  means  that 
ft,*  =  h.*,  therefore  all  <„  must  be  identically  1.  This  condition 
is  obviously  fulfilled  by  the  eigenvalues  (4.8)'. 

b)  D  *  R  (compensation  "level"  coincides  with  the  center  of  the  sphere): 
in  this  case  the  isostatic  compensation  degenerates,  such  that  the 
mean  height  becomes  independent  on  point  position  and,  therefore, 

a  constant.  As  a  consequence  only  the  zero  degree  eigenvalue  has 
to  be  equal  to  one  (it  passes  tht  operator  undisturbed)  ,  all  other 
eigenvalues  must  be  zero  (annihilation  of  all  frequencies  higher 
than  zero).  This  condition  is  also  fulfilled  by  the  eigenvalues 
(4.8)'. 

In  other  words,  (a)  represents  an  extreme  case  of  a  high-pass  filter,  (b) 
an  eytreme  case  of  a  !ow-pass  filter.  Our  isostatic  model  as  well  as  the 
actual  isostatic  compensation  will  be  located  somewhere  in  between  (a) 
ana  (b).  The  following  graph  presents  the  behavior  of  the  eigenvalues 
for  tb  ee  generally  discussed  and  used  compensation  depths,  D  =  24, 


32,  and  40  km. 
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Let  now  the  topography  he  given  in  terms  of  a  series  of  normalized 
harmonics  , 

h(')  *  £  h«  ♦-  ('>•  l4-9a> 

h*m 

then  the  corresponding  smoothed  topography  is  obtained  through 

ft(0  ■  X  <-h—  (4.9b) 

and  the  residual  topography  Ah  ■  h  -  ft 

Ah(p)  *  Y  (l-<Jh.w  (4.9c) 

<cs? 

The  corresponding  autocovariance  functions  of  h,  ft,  and  Ah  are  given  by 

cov(h,h)  a  X  h„  P„  (cosif>)  (4.10a) 

■  >0.« 

cov(ft.ft)  »  X  *2  h„P.  Ccosip)  (4.10b) 

»>0.« 

cov(Ah.Ah)  *  £  (1-k.)2 h.P.  (cos4>)  (4.10c) 

*>Q,m 

with  the  degree  variances  h.  =  S  h2 w-l . 

The  behavior  of  (1  2  for  three  compensation  depths  is  shown  in  Fig. 

4.2.  The  energy  in  the  low  frequency  part  is  dampened  because  the  low 
frequency  content  of  h  and  n  is  almost  the  same.  The  energy  in  the  high 

r\j 

frequencies,  however,  is  hardly  reduced  since  n  has  hardly  any  power  in 
the  high  frequencies.  The  deeper  the  compensation  level,  the  more  energy 
remains  from  the  high  frequent  part. 
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Fig.  4.2  Energy  dampening  factor  for  residual  height. 

What  we  are  primarily  concerned  about,  is  the  b  =  const.  -  problem. 
Considering  again  the  ideal  case  of  homogeneous  density  and  terrain- 
corrected  free-air  anomalies,  we  conclude  from  the  model  (4.1)  that  b 
is  constant  if  the  mean  height,  as  defined  by  (4  2),  is  used.  However, 
In  practical  determinations  of  b,  valid  for  a  specific  area,  a  constant 
mean  height  is  used.  In  gneral ,  ft  as  defined  by  (4.2)  is  not  constant 
over  a  limited  area  such  as  1°  x  1°.  The  error  which  we  commit  by  us¬ 
ing  a  constant  mean  height  instead  of  a  variable  mean  height  h  enters 
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fully  in  the  determination  of  b.  Denoting  the  constant  height  of  a  sped 
fic  region  by  F  and  2ttGp  by  bo,  the  factor  b,  as  determined  by  a  least  - 
squares  adjustment  (equ.  (3. 8}'’),  is  obtained  by 


b 


h,T  h' 


(4.11) 


where  h'  and  Agf  denote  the  vectors  of  residual  elevation  and  (terrain 
corrected)  free-air  anomalies, 


h,T  :  *  (hi-h,h2-F,  ...,h„-F), 

Ag':  =  (Agi-Ag,Ag2-Ag,  ....  Ag.-Ag). 


According  to  our  assumption,  the  mean  anomaly  Ag  is  given  by 

Ag  3  bo  (h  -  ft)  (4.12) 


and  therefore,  the  reduced  anomalies  can  be  represented  by 
Agr  =  b0(h  -  ft)  -  bo(h  -  ft) 

3  b0(h  -  FI  -  bo(ft  -  ft) 

■  bo  n  -  boh  .  (4.13) 


/*S 

Introducing  this  relation  ''nto  (4.11),  the  error  5b  :  3  b  -  b0  is  given 
by 


5o  =  bo 


h 

h 


'Tft' 


rT 


hf 


and,  v.'th  the  triangle  inequality  |a  bj<Ja||b|,  we  obtain  an  estimate 


5b  <.  bo 


(4.14) 
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From  (4.13)  It  is  obvious  that  the  error  5b  vanishes  If  the  size  of  the 
area,  for  which  F  and  £g  are  constant,  goes  to  zero.  Vice  versa,  6b  will 
increase  with  increasing  area  size.  Equation  (4.14)  shows  that  fib  is 
proportional  to  the  rms  variation  of  the  residual  mean  elevation  as 
defined  by  (4.2),  and  indirect  proportional  to  the  ^ms  variation  of  the 
residual  elevation.  Therefore,  an  estimate  of  b  for  small  blocks  in 
mountainous  areas  should  give  small  estimation  errors,  provided  our  model 
is  correct  and  the  data  distribution  is  sufficiently  homogeneous  and 
dense.  Poor  estimates  have  to  be  expected  for  large  blocks  in  flat  areas. 
Both  phenomena  have  been  strongly  confirmed  by  practical  determinations 
of  b  (Uotila,  1967a, b;  SUnkel  and  Malits,  1981). 

The  mean  elevation  ft  is  defined  as  a  weighted  average  of  the  ele¬ 
vation  h  with  distance-dependent  weights;  ft(r)  is  a  smooth  surface. 

This  smooth  surface,  however,  is  approximated  by  a  step  function  F  in 
all  practical  applications.  As  we  have  seen  above,  the  variation  of  ft 
with  respect  to  a  constant  is  responsible  for  the  error  in  fib,  provided 
our  model  is  valid;  the  block  size  is  closely  related  to  the  error  fib. 
Since  we  shall  hardly  ever  work  with  a  mean  surface  ft  but  rather  with  F 
we  are  interested  in  the  effect  of  fib,  caused  by  the  replacement  of  ft 
by  h.  We  will  again  consider  Faye  -  anomalies  (terrain  corrected  free- 
air  anomalies)  which  are  linearly  related  to  the  elevation  h, 

Ag  -  b0(h  -  ft)  ; 

if  we  replace  h  by  F,  we  make  b  variable  (b0  is  a  global  constant)  and 
dependent  on  position, 

Ag  =  b  (h  -  F)  . 
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The  right  hand  sides  of  both  equations  havo  to  be  equal  for  all  points, 
and  we  obtain  the  condition 

bo Ch  -  ft)  *  b(h  -  TH ,  (4.15) 

Splitting  b  up  into  b  *  b0  -  £b  and  adding  zero  to  the  right  hand  side, 
we  obtain 

bo (h  -  ft)  ■  (bo  -  <Sb)(h  -  ft  +  ft  -  IT), 
and  consequently 

6b(h  -  F)  =  b0(ft  -  F).  (4.15)' 

(Note  that  5b,  h,  %  are  variable,  F  is  constant.)  If  b  is  determined 
ny  means  of  least-squares  adjustment  with  parameters,  the  function 
5b(h  -  F)  is  considered  as  noise;  in  our  model  this  noise  is  represented 
by  the  difference  between  the  mean  elevation  surface  ft  and  the  constant 
F.  Note  that  in  the  least-squares  collocation  solution  bo (ft  -  F)  is 
treated  as  a  signal,  but  not  as  noise.  As  a  matter  of  fact,  the  power 
of  this  noise  is  a  measure  for  the  error  of  estimation  of  b  in  the  least- 
squares  adjustment  concept. 

In  the  sequel  we  shall  investigate  the  average  deviation  of  F  from 
ft.  The  derivations  are  particularly  simplified  if  F  is  considered  as 
the  mean  elevation  o.er  a  circular  region  and,  moreover,  if  F  is  consid¬ 
ered  as  the  output  of  a  moving  average  applied  to  the  actual  topography 
h.  These  two  simplifying  assumptions  do  not  quite  reflect  reality,  but 
the  advantage  of  working  with  the  concept  of  isotropy  justifies  this 
fr.mal  inaccuracy. 
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The  moving  average  IT  of  h,  taken  over  a  circular  cap  with  radius 


,  can  be  expressed  by 


^  8nhnm<j)nm(>) 
n,m 


(4.16) 


for  the  same  reason  and  in  the  same  way  as  h  was  expressed  by  (4.9b). 
In  this  context  Bn  ar®  the  eigenvalues  of  an  isotropic  moving  average 
operator  with  an  integral  kernel  B(t;  to)  defined  by 


,  Ml -to)  for  toitil 

0  else 

(t  denotes  the  cosine  of  the  spherical  distance.)  The  eigenvalues 


(4.17) 


are  obtained  through 


en(t0) 


1-to  J  dt* 


The  following  expression  can  be  found  in  (Meissl,  1971,  p.  24): 
6n(to)  3  Pt7  IthT  [  Pn-iCto)  '  Pn+i(ta)]  ; 


(4.18a) 


Sjoberg  (1980)  has  recently  derived  a  quite  attractive  recurrence  rela¬ 
tion  which  does  not  require  the  computation  of  the  Legendre  polynomials: 

Bn(to)  =  —j  [(2n-l)to6n^(t0)-  (n-2)6n_2(t0)  ],  na>2 


So  =  1  and  Bi(t0)  »  H  (1+to). 


(4.18b) 


With  (4.9b)  and  (4.16),  the  relative  error  5b/b0  (equ.  (4.15)')  can 


be  represented  by 
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»L.f, \  ^  (Kn”®n)hnm^nm(f^ 

Mil  3  Hi™ _ 

*)°  X  (l-Sn^nm^nm^^ 
n,m 


(4.15)" 


It  is  obvious  that  6b  vanishes  if  Kn=3n  for  all  n.  Due  to  the  different 
behavior  of  Kn  and  $n,  this  condition  is  only  fulfilled  in  the  extreme 
cases  of 

a)  D  =  0  and  ipo  s  0  , 

b)  D  =  R  and  \pc  -  ^  . 

For  all  other  realistic  situations  like  D=24,  32,  or  40  km  as  discussed 
here,  Bn  will  in  general  be  different  from  therefore,  b  will  not 
vanish.  The  variance  of 


S(Kn"^n^nm  ^nm^^’ 
n,m 

given  by  2(Kn-Bn)2hn  (4-19) 

n 


is  obviously  a  measure  for  the  mean  square  value  of  6b.  Therefore,  the 
goal  is  to  minimize  (4.19)  which  can  be  achieved  by  selecting,  for  each 
compensation  depths  0,  an  appropriate  cap  radius  tp0-  However,  the 
optimal  relation  between  D  and  ip0  is  influenced  by  the  actual  degree 
variances  of  the  topography.  This  means  that  we  have  to  know  the  degree 
variances  hn  of  the  topography.  Rapp  (1981)  has  recently  determined  the 
degree  variances  of  rock  topography  for  n£l80  based  on  a  1°  x  1°  mean 
elevation  data  set.  They  can  be  used  in  (4.16)  to  represent  the  energy 
which  is  contained  in  the  long-to-medium  wavelengths  of  the  topography. 
Higher  degree  variances  have  to  be  obtained  from  a  degree  variance  model. 
Since  the  author  is  not  aware  of  the  existence  of  an  appropriate  model. 


37 


another  way  has  been  chosen:  In  the  medium  frequency  range  the  actually 
"observed"  free-air  gravity  anomaly  degree  variances  and  the  correspond¬ 
ing  ones  derived  from  topographic  data  agree  fairly  well;  this  agreement 
should  be  even  better  in  the  high  frequency  range  because  of  the  strong 
linear  correlation  between  gravity  and  elevation.  T;ierefore,  it  was 
quite  natural  to  use  this  source  of  information  as  topographical  height 
degree  variance  model  for  medium  to  high  frequencies. 

Two  gravity  anomaly  degree  variance  models  which  fit  real  world 
gravity  data  best  (potential  coefficients  to  degree  180  based  on  a  com¬ 
plete  set  of  1°  x  1°  mean  free-air  anomalies,  an  observed  variance  of 
1800  mgal2,  and  a  variance  of  the  horizontal  gravity  gradient  of  800  E2), 
are  two  parameter  models  suggested  by  Moritz  (1977),  numerically  investi¬ 
gated  by  Jekeli  (1978)  and  Rapp  (1979).  Both  models  have  the  form 


Cn  =  (n-1) 


(4.20) 


and  are  determined  by  6  parameters  each.  '‘Case  One"  model  of  Rapp  (1979) 
gives  the  best  overall  fit  to  the  data,  "Case  Two"  model  fits  the  observed 
degree  variances  best. 


model 

ai£mgal 2] 

c(2  [mgal 2] 

Ol 

a2 

Ai 

a2 

"Case  One" 

3.4050 

140.03 

.998006 

.914232 

1 

2 

"Case  Two" 

14.966 

999.25 

.987969 

.850000 

75 

20 

Table  4.1  Used  "best"  degree  variance  model  parameters, 
from  Rapp  (1979) ,  p.  15. 


(meter) 
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For  a  linear  correlation  model  between  gravity  and  t  levation  with  the 
Bouguer-factor  2irGp  =  0,112  mgal/m,  the  relation  beuren  the  degree 
variances  hn  (elevation)  and  cn  (gravity)  is  given  by 


Figures  4.3a,b  show  the  degree  variances  of  observed  rock  topography  up 
to  degree  180,  and  such  ones  derived  from  Rapp's  "Case  One  &  Two"  gravity 
anomaly  degree  variance  models  for  n>180. 


4 

0 


degree  n 

Fig,  4.3a  Observed  (n i  180)  and  from  (4.20)  and  (4.21) 

derived  (n>180)  rock  topography  degree  variances 
("Case  One"  -  model)  based  on  various  compensa¬ 
tion  depths. 
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degree  n 

Fig.  4.3b  Observed  ( n  <.  180)  and  from  (4.20)  and  (4.21) 

derived  (n>180)  rock  topography  degree  variance^ 
("Case  Two"-model)  based  on  various  compensation 
depths. 


It  is  obvious  from  the  above  two  figures  that  the  "Case  Two"-mode1  has 
much  less  power  in  the  high  frequency  part  than  the  "Case  One"-model. 
"Case  One"-model  matches  the  trend  obviously  much  better;  the  50  km 
compensation  depths  seems  to  be  optimal  for  the  "Case  One"  model,  32 
km  for  the  "Case  Two"-model. 
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The  contribution  of  each  degree  variance  to  (4.19)  is  weighted  by 
(<n-8n)2.  These  weights  are  graphically  represented  in  Figures  4.4a-c 
for  various  compensation  depths  and  cap  sizes.  Needless  to  say,  the 
zero  line  corresponds  to  the  never  toltilied  ideal  <_ase  =  6n,  vn. 

The  following  conclusions  can  be  drawn  from  these  figures:  The  weights 
depend  strongly  on  the  assumed  compensation  depth;  for  the  cases  con¬ 
sidered  here  (0=24,32,  and  40  km),  the  cape  sizes  tj>o=15'  and  ^0a2o  can 
be  ruled  out  because  of  their  too  large  deviation  from  the  ideal  case; 
for  the  generally  adopted  compensation  depth  Da32  km  the  overall  weight 
minimum  is  somewhere  around  tp0=30'  to  60'.  The  weights  are  significantly 
different  from  zero,  up  to  degree  n=750  which  corresponds  to  a  wavelength 
of  about  50  km;  high-frequent  variations  (n>750)  of  the  topography  are 
obviously  very  similarly  averaged  by  (4.2)  and  (4.17);  very  low-frequent 
variations  (n<36)  are  similarly  averaged  by  (4.2)  and  (4.17)  alike;  the 
difference  between  Kn  and  6n  is  strongest  in  the  medium- frequency  range. 
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Since  (4,16)  does  not  only  depend  on  <n  and  3n  but  also  on  hn,  it  is 
particularly  important  to  have  a  good  estimate  for  medium-degree  rock 
topography  degree  variances  available.  Observed  rock  topography  degree 
variances  up  to  degree  180  were  available  to  the  author  due  to  the 
excellent  correspondence  between  observed  and  anomaV-derived  degree 
variances  for  n  =  120  to  180,  particularly  for  the  "Case  One"-model  of 
Rapp  (1979),  it  was  decided  to  use  this  model  as  a  representative  one 
for  dfc.  <*es  n>180.  With  this  data  we  computed  first  the  L2-norm  of 

bo||ft-F||  *  | <6b,  h-h>|  (4.22) 

with  II  f  H2  :=  da  . 

a 

for  various  compensation  depths  0  and  cap  radii  <P0.  Naturally,  the  norm 
depends  on  D  as  well  as  on  fy,;  for  the  optimal  choice  of  ,  corresponding 
to  a  prescribed  D,  the  norm  depends  strongly  on  i^0,  but  weakly  on  D,  and 
assumes  values  between  8  and  11  mgal;  (this  corresponds  to  a  minimum  rms 
difference  between  ft  and  Ti  of  some  70  to  100  meters.)  The  8  mgal  value 
has  been  obtained  for  D=24  km,  the  "Case  Two"-  model  and  a  cap  radius 
of  tj;0=45'.  Due  to  Schwarz'  inequality  we  are  able  to  estimate  a  globally 
valid  lower  bound  of  the  error  6b  through 

||ft-h  || 

|| 6b  ||  ^  b0  — —  .  (4.23) 

II  h-h  || 

For  the  optimal  choice  of  '^0,  corresponding  to  a  prescribed  compensation 
depth  0,  the  lower  bound  of  ||5  b  ||  is  practically  constant  and  equals 
3.0*10’2  which  is  30%  of  the  normal  Bouguer  factor.  Those  optimal  values 
have  been  obtained  for  0=24  km,  ^0=40'  down  to  0=40  km,  ^0=110"  and  the 


"Case  One"-  model.  The  graphs  fn  Figs.  4.5  and  4.6  show  b0||  h-K  ||  and 
||  5b  ||  as  defined  by  (4.23}  dependent  on  various  compensation  depths 
and  cap  sizes,  for  both  the  "Case  One"  and  the  "Case  T*o"  -  model. 


•=  10  1 


Case  One  -  model  (n  >  180) 


bollh  -  h  ll/ll  h  -F  II  (mgal/m) 
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In  the  previous  investigation  of  the  error  in  b,  caused  by  a  replacement 
of  h  by  h,  only  the  global  case  has  been  considered.  Therefore,  the 
optimal  estimates  for  t|>a  should  be  interpreted  with  this  reservation. 

Local  optimal  ip0  -  values  can  be  considerably  smaller;  therefore,  the 
<|>o  -  estimates  represent  rather  upper  limits.  Best  local  estimates  de¬ 
pend  on  the  Individual  situation;  they  could  be  obtained  on  the  basis 
of  a  detailed  digital  terrain  model. 

Summarizing  we  can  say  that  large  deviations  of  b  from  the  normal 
Bouguer  gradient  of  0.112  mgal/m  can  be  expected  if  b  is  determined  for 
a  large  block  by  least-squares  adjustment  and  in  addition,  if  the  Bouguer 
anomaly  is  not  constant  within  the  block.  Using  a  simplified  concept 
of  isostatic  compensation,  we  could  give  a  very  simple  mathematical  ex¬ 
planation  for  this  phenomena.  If  the  least-squares  adjustment  concept 
is  used  for  the  estimation  of  b,  the  selection  of  an  appropriate  block 
size  should  be  done  very  carefully.  In  general  it  is  much  better  to 
choose  a  too  small  block  size  than  a  too  largo  one;  this  is  particularly 
true  for  mountainous  areas.  The  least-squares  collocation  determination 
of  b  is  quite  insensitive  with  respect  to  the  choice  of  the  block  size 
because  it  takes  into  account  the  variation  of  the  Bouguer  -  anomaly 
within  the  block.  Therefore,  it  is  a  very  good  advise  to  estimate  b  using 
the  method  of  least-squares  collocation  with  parameters,  particularly 
in  areas  with  sparse  data  coverage.  In  addition,  collocation  allows  at 
the  same  time  the  estimation  of  the  Bouguer  anomaly  field  and  even  the 
prediction  of  surface  free-air  point  and  mean  anomalies  together  with 
their  accuracies.  The  author  is  not  aware  of  any  other  nearly  as 
powerful  existing  method. 
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5.  THE  TERRAIN  EFFECT  ON  POINT  ANOMALIES 

The  natural  goal  of  estimation  problems  is  to  keep  the  error  of 
estimation  as  small  as  possible.  Least-squares  interpolation,  in  parti¬ 
cular,  is  sensitive  with  respect  to  the  statistical  properties  of  the 
field  to  be  estimated,  expressed  in  terms  of  a  (usually)  isotropic  and 
homogeneous  covariance  function.  The  interpolation  error  depends  strongly 
on  a)  the  variance  Co  and  b)  on  the  ratio  r  «  correlation  length/data 
spacing  (Siinkel,  1981).  A  small  interpolation  error  is  achieved  by  a 
small  Co  and  a  large  r.  Therefore,  any  data  reduction  process,  which 
decreases  the  variance  and  increases  the  correlation  length,  has  to  be 
favorably  considered.  It  is  cannon  sense  that  the  irregularities  of  the 
topography  account  significantly  in  the  power  of  first  and  higher  order 
derivatives  of  the  gravitational  field.  In  particular,  the  free-air 
anomaly's  high-frequent  variation  comes  primarily  from  the  influence  of 
the  topography.  In  other  words,  the  topography  makes  C0  increase  and 
the  correlation  length  decrease,  fhis  is  why  predictions  in  mountainous 
areas,  based  on  unreduced  free-air  gravity  field  quantities,  give  poor 
accuracies.  If  we  want  to  achieve  high  Drediction  accuracy,  we  have 
only  two  alternatives:  do  manpower-consuming  expensive  fio1d  work  and 
collect  more  data  just  to  make  r  increase,  or  reduce  the  data  for  the 
influence  of  the  topography.  Needless  to  say,  plain  mortal  geodesists 
prefer  the  latter. 

In  linear  and  planar  approximation  the  topographic  correction  of 
gravvty  at  a  surface  point  P  is  given  by  (Moritz,  1969,  p.  10) 
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C,  -  h  GpR^l^M'dq  (5.1) 

o 

with  Ho  *  2R  sin  2.  Since  the  fntegral  kernel  Ho3  drops  rapidly  to 
zero,  only  a  very  small  region,  centered  at  the  computation  point  P,  has 
to  be  considered  for  the  evaluation  of  (5.1);  therefore,  it  is  legiti¬ 
mate  to  formally  replace  the  sphere  by  its  tangential  plane  at  P. 

In  the  following  we  shall  investigate  how  detailed  topographic 
information  has  to  be  made  available  in  order  to  meet  certain  accuracy 
requirements.  Let  us  first  investigate  the  critical  zone  in  the  neighbor¬ 
hood  of  the  computation  point  P.  following  Heiskanen  &  Moritz  (1967, 
p,1214i$.)»  we  represent  the  topography  around  P  in  terms  of  a  Taylor 
series, 


h(s,ot)  «  h*  +  s(hx  cosa  +  hy  sina)  +  ...;  (5.2) 

here  s  and  a  denote  the  planar  distance  and  the  azimuth,  x  and  y  are 
cartesian  coordinates.  Then  (5.1)  is  represented  by 

27T  So 

5Cp  =  *$  Gp  j  J s  ds  da  , 
a=o  Saa 


and  with  (5.2)  we  obtain,  neglecting  higher  order  terms. 


6Cp  a  H  Gp 


// 


(h£  cosa+  hy  sina+  2hxhy  sina  cosa)dsda. 


a=o  s= o 

Due  to  the  orthogonality  relations  of  trigonometric  functions,  this  ex 
pression  reduces  to  the  simple  form 
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<5Cp  *•  trSps0»s(h2x  +  hy).  (5.3) 

Consider  e.g.  a  small  zone  with  a  radius  s.  =  300  m  and  a  moderate  slope 
of  20°  only,  the  terrain  correction  will  assume  a  value  of  about  1  mgal . 
Consequently,  the  resolution  of  the  used  terrain  model  has  to  be  very 
high  in  the  neighborhood  of  the  calculation  point,  unless  P  is  located 
on  a  "flat  spot"  of  the  terrain. 

In  order  to  study  the  response  of  C  to  the  terrain,  we  need  terrain 
models.  However,  the  choice  of  a  proper  terrain  model  is  a  very  delicate 
problem.  To  some  degree  it  can  be  anticipated  that  C  will  be  relatively 
insensitive  with  respect  to  high-frequent,  and  sensitive  with  respect  to 
medium-frequent  topographic  variations.  As  a  matter  of  fact,  C  should 
depend  somehow  on  the  power  of  those  variations.  For  the  sake  of  sim¬ 
plicity  we  choose  a  very  simple  but  instructive  topographic  model:  an 
isotropic  model  represented  by 

h(s)  *  h0  cos  w  s,  (5.4) 

centered  at  th<*  computation  point  P;  (note  that  hP  =  hQ. )  The  model 
looks  quite  unnatural,  but  it  isn't:  imagine  a  gravity  station  either 
on  a  top  of  a  mountain  surrounded  by  (circular)  mountain  chains  of  com¬ 
parable  height  which  are  separated  by  valleys,  or  -  numerically  equivalent- 
a  gravity  station  in  a  valley  surrounded  by  mountain  chains  and  other 
valleys.  The  only  unlikely  structure  in  our  model  is  the  circular  sym¬ 
metry,  constant  amplitude  and  frequency.  However,  in  order  to  study  the 
influx-  o  of  topography  on  gravity,  sacrifices  have  to  be  made  resulting 
in  simple  structures. 


Wtth  (.5,4)  the  topographic  correction  (5.1)  is  given  by 


(5.5) 


(The  integration  over  the  azimuth  has  already  been  performed.)  According 
tn  Ryshik  and  Gradstein  (1963,  p.  114,  No.  2.523;  p.  115,  No.  2.526)  the 
integral  in  (5.5)  assumes  the  form 


Sso 


The  second  expression  [•]  can  easily  be  shown  to  vanish  by  a  Taylor  ser¬ 
ies  evaluation  at  s=0.  Considering 
00 

/sin  px  _  jr_ 
px  “  2p 

o 

(Ryshik  and  Gradstein,  1963,  p.  169,  No.  3.522),  the  first  expression 
(•)  assumes  the  simple  form  tt/4u>;  with  u>  =  2irA  wave-length), 

the  topographic  correction  (5.5)  is  given  by 
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This  remarkably  simple  result  for  the  equally  simple  model  deserves  a 
discussion:  Let  us  first  reflect  about  its  validity.  Formula  (5.1) 
represents  the  linear  term  of  a  series  expansion  of  the  topographic 
correction  in  planar  approximat  or  with  respect  to  [(h-h, )/£<>] 2 »  which 
essentially  represents  the  square  of  the  tangent  of  the  elevation  angle 
3  of  a  variable  terrain  point  with  respect  to  P,  tan2g.  The  series 
obviously  converges  if  |8j<  45°;  therefore,  equations  (5.1)  as  well  as 
(5.5)' are  valid  approximations  for  moderate  terrain  only  with  tan26 

i 

significantly  smaller  than  1.  In  terms  of  the  wavelength  A  of  our 
model,  this  translates  into  A>  4h0;  in  order  to  be  absolutely  save,  we 
should  rather  write 

h2 

C  =  *3Gp  ^  ,  X  » 4h0  .  (5.5)" 

C  vanishes  if  A-*«  as  it  should  be.  The  most  important  result  (at 
least  for  the  model  considered  here)  can  be  summarized  as  follows: 

The  topographic  correction  is  proportional  to  the  square  of  the 
amplitude  h0  (and  consequently  linearly  dependent  on  the  variance  of 
the  topography),  u  J  inverse  proportional  to  the  wavelength  A  of  the 
topography,  provided  A  »4h0. 

The  following  Table  5.1  shows  C  for  our  model,  dependent  on  various 


choices  of  ho  and  A. 
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\X(m) 

100 

1000 

5000 

10  000 

10 

0.55 

0.05 

0.01 

<0.01 

100 

5.5 

1.1 

0.55 

500 

27.6 

13.8 

1000 

55.2 

Table  5.1  Topographic  correction  (mgal)  dependent 
on  amplitude  h0  and  wavelength  X. 


The  figures  obtained  with  our  simple  mod°1  agree  remarkably  well  with  real 
world  observations;  cf.  Heiskanen  and  Vening  Meinesz  (1958,  p.  154).  (Note, 
for  example,  that  the  figure  in  the  last  line  and  last  column  corresponds 
to  the  situation  of  a  mountain  top  2000  m  above  the  surrounding  valley  and 
having  a  circular  basis  of  20  km  in  diameter.)  It  is  quite  instructive  to 
compare  these  values  with  figures  derived  from  a  cone  model  with  height  = 

2  h0  and  base-radius  =  X/2;  the  topographic  correction  can  easily  be  shown 
to  equal  C  -  8TrGpho/X;  therefore,  our  two  models  differ  only  by  8/tt2  =  0.8 
and  consequently,  the  cone-model  values  corresponding  to  the  figures  in  Table 
5.1  are  only  20%  smaller;  a  very  astonishing  result,  which  both  convinces 
the  above  rule  of  thumb  and  indicates  the  moderate  influence  of  remote  zones. 

Let  us  still  investigate  another  isotropic  model  represented  by 


h(s)  =  hQsinu)  s. 


(5.6) 


centered  at  the  computation  point  P;  (note  that  hP  =  0.)  With  (5.1)  the 
topographic  correction  is  given  by 


C  =  irGpho 


f  sin2  w 

i 


ds 
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with  (Ryshik  and  Gradstein,  1963,  p,  166,  No.  3.512)  and  w  s  2ir/X  as  before, 
we  obtain 

C  *  tt3Gp  M  ^  x>4  hfi  (5.7) ' 

which  is  identical  to  (5.5)~. 

As  a  last  model  we  will  consider  a  non- isotropic  so-called  "two- 
dimensional"  model  which  has  a  constant  profile  in  one  direction, 

h(x)  »  h0sin<ux;  (5.8) 

the  computation  point  is  supposed  to  be  located  at  x  *  0.  W i th  (5.1)  the 
topographic  correction  is  given  by 


(5.9) 


Due  to  the  symmetry  of  the  integrand  with  respect  to  x  -  0  and  y  3  o,  C  can 
also  be  expressed  by 


00  00 


x*0  y*0 


The  integration  with  respect  to  >  is  very  simple  and  yields 


and  therefore  C  reduces  to 


oo 


55 


,  .  ,r.k2  /  sin2  ux 

C  -  3  Gph0  J  —jr-  dx  , 

x=Q 

which  differs  from  (5,7)  only  by  a  factor  2/tt  and  we  obtain  with  to  =  2tt/x 
as  above 

h2 

C  =  2ir2Gp  ^  ,  X»4  h0.  (5.9)" 

Also  for  this  quite  different  model  (note  that  the  computation  point  is 
located  at  the  slope  of  a  mountain  chain)  we  observe  the  same  dependence  on 
ho  and  X.  Therefore,  we  conclude  that  our  rule  of  thumb 

C  =  0 ( h o / X )  (5.10) 

seems  to  be  of  general  validity. 

What  is  the  logical  consequence  for  the  design  of  a  digital  terrain 
model  for  the  purpose  of  reducing  gravity  measurements?  Observing  (5.10), 
the  obvious  answer  is  as  follows:  the  terrain  sampling  rate  must  be  selected 
terrain-specifical ly;  rough  terrain  requires  higher  rates,  smooth  terrain 
lower  sampling  rates;  it  is  not  a  good  advice  to  keep  the  sampling  rate 
globally  constant.  If  we  are  talking  about  a  rough  terrain  we  mean  not  only 
small  wavelengths  but  also  high  power;  high-frequent  terrain  oscillations 
with  a  small  amplitude  need  not  be  resolved,  they  should  be  smoothed  out  by 
an  appropriate  smoothing  process;  a  mean  value  representation  over  rectangular 
blocks  or  even  a  mean  value  reproducing  smooth  representation  should  be 
favorably  considered.  The  block  size,  in  turn,  depends  strongly  on  the 
power  of  the  high-frequent  portion  of  the  terrain  spectrum,  and  in  addition 
on  the  desired  accuracy  of  the  terrain  correction. 


c'si^f  iH&sin4,d*-  (5*ur 

1^=0 

Observing  the  rapid  decrease  of  che  function  sin^MoW.  expression  (5.11)' 
can  be  considerably  simplified  for  practical  applications;  with 

2M-  =  2  sin  tp/2  -cos^/2  •  1 
loUT  8  sin  p 


iifrEAin-nf  j;  i 
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Oi>0  a  1°5  is  sufficient  in  any  case  (Hay-Ford  zone  0);  cf.  (Heiskanen 
and  Vening  Meinesz,  1358,  p.  154),)  In  terms  of  the  linear  distance 
s  3  fty,  above  formula  is  written 

so 

C  =  ttGo 

s=0 


/ 


ds. 


(5.11) 


Let  us  now  recall  the  isotropic  terrain  model  (5.4)  discussed  earlier. 
Comparing  (5.5)  with  (5.11)"  ,  we  observe  that  our  requirements  can 
be  considerably  relaxed:  the  model  has  to  be  such  that  only  its  mean 
square  value  (extended  over  the  azimuth  a;  s  fixed)  behaves  like 
ho  (cos  ^  s-1 ) 2 ,  otherwise  it  is  largely  arbitrary.  Consider  the  in- 

A 

tegrand  of  (5.5):  with  a  constant  denominator  the  maximum  would  be 
attained  at  s  =  (2k+l)X/2,  k  *  0,  1 .  Due  to  the  rapidly  increas¬ 

ing  denominator  s2,  the  maximum  of  the  composite  integrand  is  shifted  to 
smaller  values  of  s;  in  addition,  C  gets  its  power  mainly  from  the  re¬ 
gion  of  the  inegrand's  first  maximum;  local  maxima  of  larger  s  contribute 
very  little  to  the  terrain  effect.  This  is  why  the  simple  cone  model 
discussed  above  differs  only  little  from  the  model  (5.4)  as  far  as  the 
terrain  effect  is  concerned.  Fig.  5.1  illustrates  the  behavior  of  the 
integrand  corresponding  to  model  (5.4)  for  various  wavelengths  (common 
h0,  normalized).  It  is  quite  remarkable  that  these  graphs,  despite  of 
the  underlying  model's  simplicity,  agree  favorably  well  with  graphs 
derived  from  real  world  topography;  cf.  (Mathisen,  1976);  this  fact 
confirms  once  more  the  h2/X  -  law. 

The  practical  consequences  are  as  follows:  the  location  of  the  first 
maximum  is  strongly  influenced  by  the  variation  of  the  topography  in  the 
neighborhood  of  the  computation  point;  therefore  a  very  detailed 


distance  s(meters) 

Fig.  5.1  Integrand  [ho(cos  ~~  s  -  l)/s]  2  for 
various  wavelengths  X. 

digital  terrain  model  (DTM)  is  required  in  this  region,  unless  the  compu 
tation  point  is  located  at  a  flat  spot  of  the  topography.  The  degree 
of  resolution  of  the  2mM,  therefore,  depends  on  three  essential  factors: 

a)  on  the  variance  of  the  terrain  surrounding  the  computation  point  P, 

b)  on  the  wavelength  of  the  terrain  surrounding  P, 

c)  on  . ;ie  location  of  the  first  maximum  cf  the  integrand  in  (5.11)"  . 
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Table  (5.1)  and  Fig.  (5.1)  provide  us  with  a  rough  guideline  for 
a  proper  choice  of  the  required  DTM's  resolution;  however,  the  decision 
must  be  made  on  individual  circumstances. 
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6.  THE  EFFECT  OF  TOPOGRAPHY  ON  MEAN  ANOMALIES 

Mean  free-air  anomalies  as  defined  by  (2.2)  can  be  split  up  into 
two  components:  the  mean  Faye  anomaly  plus  a  mean  terrain  correction. 
Chapters  2  and  3  deal  with  the  estimation  of  point  and  mean  Faye  anoma¬ 
lies;  chapter  5  considers  essential  aspects  concerning  the  calculation 
of  point  terrain  corrections.  As  a  matter  of  fact,  the  effect  of  topo¬ 
graphy  on  mean  anomalies  equals  the  mean  terrain  effect;  its  estimation 
will  be  discussed  in  the  sequel. 

The  most  straightforward  way  of  estimating  the  mean  terrain  effect 
is  to  calculate  a  dense  grid  of  point  terrain  effects  and  take  the 
.average.  The  calculation  of  the  point  terrain  effect  requires  the 
evaluation  of  an  integral  formula  like  (5.1),  an  expensive  task  even  for 
high-speed  computer.  Note  that  equation  (5.1)  is  of  "differentiation 
type"  which  can  be  clearly  seen  by  observing  its  behavior  in  the  inner¬ 
most  zone  surrounding  the  computation  point  P.  The  mean  terrain  effect 
is  obtained  by  integrating  the  point  terrain  effects.  Do  we  really  have 
to  suffer  from  the  instabilities  of  differentiation  first  before  we  can 
enjoy  the  stability  of  integration?  There  must  be  a  direct  and  inexpen¬ 
sive  solution  to  the  problem  of  mean  .rain  correction  estimation  which 
avoids  the  up's  and  down  s  or'  differentiation  +  integration.  (A  very 
similar  problem  is  encountered  in  connection  with  the  astrogeodetic 
determination  of  the  geoid:  plumb-line  curvature  correction  and  ortho¬ 
metric  reduction;  see  Heiskanen  &  Moritz,  1967,  pp.  20C,  201.) 

!oi  us  consider  the  terrain  correction  as  given  by  equation  (5.11)'. 
The  mean  terrain  correction  C  =  M{C}  is  obviously  obtained  by  averaging 
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the  point  terrain  corrections  within  the  area  A a  of  consideration, 


C  = 


(6.1) 


Let  us  try  to  formulate  the  procedure  of  calculating  C: 

a)  keep  a  circle  with  radius  “infixed; 

b)  keep  a  point  P  e  A  0  fixed  and  make  it  the  center  of  the  circle; 

c)  calculate  the  mean  square  height  differences  as  defined  by  (5.12) 
for  this  particular  circle; 

d)  change  P  and  repeat  (b)  and  (c)  until  Aa  is  covered; 

e)  calculate  the  average  of  all  outcomes  of  (c); 

f)  change  ijj  and  repeat  (a)  -  (e)  until  [O,ii»0]is  covered  with  an  appro¬ 
priate  (e.g.  1°5) ; 

g)  perform  (5.11)  . 

In  terms  of  a  formula  it  can  be  written  as 


'J'o  /  2ir  \ 

C  =  C  J-L  ffL  f  Ahadai^ 

L  R  /  \*0JJ2nJ  )  * 

ijj=o  Ac  u»o 


(6.2) 


The  expression  {•>  represents  (c)  -  (e)  and  equals  the  variance  of  height 

differences  between  two  surface  points,  separated  by  a  distance  for  the 

(1) 

individual  area  A a  in  consideration.  This  variance  is  related  to  the 
elevation  auto-covariance  function  H( ^)  as  follows: 


2TT 


1 

2ttA  a 


If 


Ahzdada  =  2  [ H0  -  H(i|>)]  ,  H0=H(0). 


a=^  A  a 


(6.3) 


(1)...  provided  A  a  is  sufficiently  larger  than  ii>0 
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Replacing  the  spherical  by  planar  expression,  we  obtain  the  very  simple 
expression 


C 


So 


S«0 


(6.4) 


In  the  following  we  will  evaluate  this  Integral  for  various  covariance 
functions,  the  generalized  Hirvonen  models  and  the  Gaussian  covariance 
function. 

a)  The  classical  Hirvonen  model 


H(s) 


(6.5a) 


where  5  denotes  the  correlation  length,  yields  an  integral  (Ryshik  and 
Gradstein,  1963,  p.  60,  No.  2.141) 


ds 

s*o 


and  we  obtain  fcr  the  mean  terrain  correction  the  value 


C  a  •n-2Gp  ~ 


(6.6a) 


h)  The  second  Hirvonen  model  has  been  discussed  in  Moritz  (1980), 
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where  the  parameter  a  is  related  to  the  correlation  length  £  through 
5  *  aVT". 


Ha -z-i its)  ds  =!k 

sz  a 


and  we  obtain  for  C  the  value 


C  *  2  VT  ttGp  . 


(6.6b) 


c)  The  last  Hirvonen  covariance  model  (Moritz,  1980)  is 


H(s)  = 


l+(f) 


ITT*  • 


(6.5c) 


where  the  parameter  a  is  related  to  the  correlation  length  5  through 
£  ■  a(2 -  l)5*.  According  to  Ryshik  and  Gradstein  (1963,  p.  82,  No. 
2.268)  the  integral  yields 


Hq  ~  H( s)  d  =  LH, 

s7  a 


and  we  obtain  for  C  the  expression 


C  =  4V2  3  -1  ttGo  ^  . 


(6.6c) 


d)  The  Gaussian  covariance  model 


H(s)  =  H0e 


-a^7- 


(6.5d) 


1 
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has  a  correlation  length  £  *  Y  v£n2  and  yields  an  integral  (Ryshik  and 
Gradstein,  1963,  p.  151,  No.  3.273) 


f  ds 

s=a 


Ho 

a 


and  we  obtain  for  the  mean  terrain  correction  the  value 


C  *  Z^TYnP  TTGp  ^  .  (6.6d) 

Lj 

In  11  four  covariance  models  we  observe  the  common  factor  irGp  . 

The  model -specific  multiplication  factor  is  shown  in  Table  6.1. 


Covariance  model 

Hirvonen  (a) 
Hirvonen  (b) 
Hirvonen  (c) 
Gauss 


Co 


IT 

=  3.14 

2  VT 

=  3.46 

4  V2'/j  -1 

=  3.07 

2  Vtt  In  2 

=  2.95 

Table  6.1  Model -speci fic  multiplication  factors, 

C  -  Co'  ttGp  ^  . 

We  notice  that  for  all  four  models  the  mean  terrain  correction  C"  varies 
in  a  very  narrow  range  with  variations  of  only  10%.  Therefore,  we  con¬ 
clude  iK.*, i  the  estimation  of  C  is  only  little  sensitive  with  respect  to 
the  model  choice.  We  summarize: 
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The  mean  tenMin  coKAection  depends  lineaAly  on  the  tevuxin  variance 
and  is  invense  paopoational  to  the.  connelation  Length  o h  the  teonain  c o- 
va/uance  function  ootid  ^on  the  comi.deA.ed  cuiea.  It  depends  weakly  on 
the  type  oh  the  cooaniance  model. 

This  very  interesting  and  astonishingly  simple  result  deserves  a 
discussion.  The  variance  of  the  height  differences  as  defined  by  (6.3) 
is  valid  for  the  area  Aa  in  consideration;  therefore,  the  terrain  covari- 
»  ance  function  H(s)  has  to  be  derived  from  topographic  data  in  the  same 
limited  area  only.  In  flat  areas  the  variance  of  the  topography  is  small 
and  the  correlation  length  big,  resuiting  in  a  very  small  mean  terrain 
effect  as  to  be  expected.  In  mountainous  areas,  in  constrast,  we  observe 
a  big  variance  and  a  small  correlation  length,  resulting  in  a  very  signi¬ 
ficant  mean  terrain  effect  as  also  to  be  expected.  How  does  it  come 
that,  according  to  equations  (6.6a-d),  the  mean  terrain  effect  is  vir¬ 
tually  independent  on  the  size  cf  the  area  Aa?  The  answer  is  simple: 
the  area  size  A a  lurks  in  the  background  and  is  implicitly  introduced 
through  the  definition  and  determination  of  the  covariance  parameters  H0 
and  5.  It  is  known  from  experience  that  the  terrain  variance  and  correla¬ 
tion  length  can  considerably  change  with  the  size  of  the  area  Aa.  If  A a 
increases,  the  correlation  length  will  in  general  increase  too.  H0  can 

attain  quite  large  values;  in  mountainous  areas  values  of  H0  =  0.5  to  1 

2 

km  can  be  observed;  the  correlation  length  of  the  terrain  covariance  func¬ 


tion  is  usually  much  smaller  than  the  one  of  the  free-air  anomaly  covari¬ 
ance  function;  particularly  in  mountainous  areas  it  can  be  as  small  as  a 
very  few  kilometers,  Fig.  6.1  shows  lines  of  equal  mean  terrain  effect 


depending  on  the  r.m.s,  elevation  and  the  correlation  length  £  for  the 
Hirvonen  (a)  -  model;  Fig.  6.2  is  a  3-D  representation  of  the  behavior 
of  (6.6a). 


correlation  length  E,  (km) 

Fig.  6.1  Lines  of  equal  mean  terrain  effect  depending  on 
the  r^.m.s.  variation  of  elevation  and  the 
correlation  length  E,  for  the  Hirvonen  (a)-model. 

unit:  mgal 
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7.  CONCLUSIONS  AND  RECOMMENDATIONS 

Different  methods  have  become  available  for  the  prediction  of  point 
and  mean  free-air  anomalies  ba>*d  on  observed  point  anomalies.  By  far 
the  most  attractive,  but  also  most  expensive,  tool  is  least-squares  pre¬ 
diction.  The  most  general  version  of  least-squares  gravity  prediction 
involves  the  processing  of  gravity  and  elevation  data  simultaneously. 

In  this  report  a  least-squares  approach  is  proposed,  which  takes  the  gen¬ 
erally  observed  strong  linear  correlation  between  (terrain  reduced)  free- 
air  anomalies  and  elevation  into  account  in  terms  of  two  model  parameters; 
the  method  is  therefore  based  on  least-squares  collocation  with  parameters. 
The  signal  is  basically  a  residual  Bouguer  anomaly  and  such  is  the  cor¬ 
responding  covariance  function.  Bouguer  anomalies  are  smooth,  therefore, 
the  covariance  function  has  a  rather  long  correlation  length;  consequently, 
the  interpolation  (prediction)  accuracy  is  high.  It  has  been  shown  that 
a  regional  variation  of  the  Bouguer  anomaly  causes  the  Bouguer  factor  to 
decrease  if  a  regional  least-square s  adjustment  is  used;  in  addition  its 
estimation  accuracy  becomes  very  poor.  The  collocation  solution,  in  con¬ 
trast,  models  the?  Pougu er  anomaly  field  statistically  and  provides  highly 
reliable  estimates  of  the  Bouguer  factor  which  turned  out  to  be  generally 
quite  close*  to  its  normal  value.  Particularly  important  for  free-air 
anomaly  predictions  in  mountainous  areas,  is  the  data  reduction  for  the 
influence  of  the  topography.  Various  simple,  but  instructive,  topograph¬ 
ical  models  have  been  studied.  It  turned  out  that  the  mean  terrain  effect 
is  lir  •  r T y  dependent  on  the  variance  and  ^'.directly  proportional  to  the 
correlation  length  of  the  residual  topography  within  the  area  of  consider¬ 
ation.  Therefore,  the  accuracy  of  its  determination  depends  on  how  good 
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the  estimates  of  the  terrain  variance  and  the  terrain  correlation  length 
are. 

Summarizing  we  propose  the  following  method  for  mean  free-air  gravity 
anomaly  estimation: 

a)  Reduce  gravity  data  for  the  terrain  effect; 

b)  determine  an  empirical  covariance  function  and  fit  a  model; 

c)  determine  the  correlation  model  parameters  (mean  Bouguer  anomaly 
and  the  Bouguer  factor)  by  least-squares  collocation; 

d)  predict  point  and/or  mean  anomalies  by  least-squares  collocation  with 
parameters ; 

e)  add  the  point  and/or  mean  terrain  effect  to  the  obtained  estimate 
in  order  to  get  free-air  anomalies. 

Predictions  with  real  world  data  have  shown  that  free-air  5'  x  5'  anoma¬ 
lies,  even  in  very  rough  mountainous  areas,  can  be  predicted  with  an  ac¬ 
curacy  of  <±5  mgal.  using  the  method  proposed  above,  provided  the  data 
density  is  better  than  1/10  km2  and  the  effect  of  topography  has  been 
carefully  taken  into  account. 

Further  investigation,  particularly  as  far  as  the  resolution  of  the 
digital  terrain  model  is  oncerned,  are  both  useful  and  necessary. 
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